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ABSTRACT 

The  3D  shape  of  a  textured  surface  is  recovered  from  its  projected  image  on  the  as¬ 
sumption  that  the  texture  is  homogeneously  distributed.  First,  the  homogeneity  of  a 
discrete  texture  consisting  of  dots  and  line  segments  is  defined  in  terms  of  the  theory  of 
distributions.  Next,  distortion  of  the  observed  texture  density  due  to  perspective  projec¬ 
tion  is  described  in  terms  of  the  first  fundamental  form,  which  is  expressed  with  respect 
to  the  image  coordinate  system.  Based  on  this  result,  the  basic  equations  to  determine 
the  surface  shape  are  derived  for  both  planar  and  curved  surfaces,  and  numerical 
schemes  are  proposed  to  solve  them.  Necessary  data  are  obtained  in  the  form  of  summa¬ 
tion  or  integration  of  functions  over  the  texture  elements  on  the  image  plane.  Ambiguity 
in  the  interpretation  of  curved  surfaces  is  also  analyzed.  Finally,  numerical  examples  for 
synthetic  data  are  presented,  and  our  method  is  compared  with  other  existing  methods. 
It  is  shown  that  all  other  methods  can  be  explained  in  terms  of  our  formulation. 
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1.  INTRODUCTION 

The  recovery  of  the  3D  shape  of  a  surface  from  its  projected  image  is  one  of  the 
most  important  challenges  in  computer  vision.  The  computation  is  based  on  various 
clues  such  as  texture,  shading  and  motion.  3D  recovery  from  texture  is  possible  if  we 
have  some  knowledge  about  the  true  texture.  If  the  projected  texture  has  different  pro¬ 
perties  from  those  of  the  true  texture,  the  3D  shape  is  recovered  from  the  difference 
between  the  observed  properties  and  the  original  properties  that  we  know.  For  example, 
if  the  true  texture  is  known  to  be  an  array  of  texture  elements  of  a  known  shape,  the 
surface  shape  can  be  inferred  from  the  observed  distortion  of  the  texture  elements. 

If  we  do  not  know  the  true  texture  but  know  its  statistical  properties,  we  can  draw 
an  inference  from  the  difference  between  the  observed  statistical  properties  and  the  true 
properties.  Assumptions  often  adopted  are  homogeneity  and  isotropy  of  the  true  texture. 
The  assumption  of  isotropy  asserts  that  the  texture  consists  of  line  segments  with  no 
preferred  orientations.  A  clue  for  3D  recovery  is  obtained  if  the  observed  texture  seg¬ 
ments  have  a  preferred  orientation.  This  approach  was  first  investigated  by  Witkin  [l], 
and  the  algorithm  was  improved  by  Davis,  et  al.  [2].  Kanatani  [3]  gave  a  rigorous 
mathematical  description  of  the  problem  and  explicit  formulae  to  solve  the  problem  by 
means  of  tensor  calculus  and  the  principle  of  stereology. 

The  assumption  of  homogeneity,  on  the  other  hand,  asserts  that  the  texture  is  uni¬ 
formly  distributed  over  the  surface.  When  projected,  the  texture  becomes  dense  on  the 
image  of  the  surface  part  away  from  the  observer  and  sparse  on  the  part  near  the 
observer.  This  clue  has  been  considered  long  since  by  people  like  Gibson  [4,  5],  Bajcsy 
and  Lieberman  [6]  and  Stevens  (7j,  but  their  argument  was  based  on  naive  intuition.  It 
was  not  until  Aloimonos  and  Swain  [8|  and  Dunn  [9]  that  the  problem  was  treated  in 
analytical  terms  based  on  the  geometry  of  perspective  projection.  However,  their  formu¬ 
lations  involve  many  unnecessary  ad  hoc  approximations  and  assumptions.  In  this 
paper,  we  show  a  mathematically  rigorous  treatment  based  on  differential  geometry  and 
the  theory  of  distributions,  taking  the  “discreteness”  of  texture  correctly  into  account. 

We  first  give  a  precise  definition  of  homogeneity  of  a  texture.  If  a  texture  consists 
of  dots  or  line  segments,  the  texture  density  is  a  singular  function  taking  the  value 
infinity  at  the  texture  dots  and  line  segments  and  0  elsewhere.  How  can  we  say  that  the 
density  is  uniform?  How  can  we  tell  that  a  given  texture  is  homogeneous?  We  will  give 
an  exact  definition  of  homogeneity  of  a  discrete  texture. 

We  next  give  an  exact  analysis  of  the  distortion  of  texture  due  to  perspective  pro¬ 
jection  in  terms  of  the  first  fundamental  form  expressed  with  respect  to  the  image  coordi¬ 
nate  system.  Our  formulation  consists  of  two  stages.  First,  we  present  the  basic  equa¬ 
tions  to  determine  the  surface  shape  for  both  planar  and  curved  surfaces  in  their  exact 
forms.  Although  they  are  difficult  to  solve  directly,  we  can  infer  various  theoretical 
consequences,  among  which  is  the  ambiguity  in  the  interpretation  of  curved  surfaces. 
We  list  all  possibilities  completely. 

Then,  we  propose  various  numerical  schemes  to  solve  these  equations,  employing 
first  order  approximation,  simulation  of  camera  rotation  and  Newton-Raphson  type 
iterations,  and  give  numerical  examples  for  synthetic  images.  Good  results  are  obtained 
even  for  a  very  sparse  texture,  and  the  estimation  approaches  the  true  value  as  the  tex¬ 
ture  density  increases. 

Lastly,  our  formulation  is  compared  with  those  of  Aloimonos  and  Swain  [8]  and 
Dunn[9].  Our  formulation  is  general  enought  to  explain  their  methods  in  our  terms. 
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Various  aspects  related  to  ■ -oplications  and  implementations  of  our  method  are  also  dis¬ 
cussed. 


2.  TEXTURE  DENSITY  AND  HOMOGENEITY 

We  consider,  in  this  section,  how  to  define  the  density  of  a  discrete  texture.  Con¬ 
sider  textures  composed  of  dots  or  line  segments.  If  we  are  to  seek  a  function  J[x,y) 
describing  the  amount  of  texture  divided  by  the  area  occupied,  we  are  forced  to  consider 
delta-function-like  singularities,  the  value  of  J{x,y)  being  infinity  at  the  texture  elements 
and  0  elsewhere,  since  the  area  of  a  dot  or  a  line  segment  is  0.  This  kind  of  singularities 
make  analytical  treatment  very  difficult.  One  way  to  avoid  singularities  is  to  regard  the 
texture  density  as  a  functional. 


DIRAC  DELTA  FUNCTION 

Let  us  recall  the  definition  of  the  (Dirac)  delta  function  6(z).  To  be  precise,  it  is  not 
a  function;  if  a  function  takes  the  value  zero  except  at  one  point,  its  integral  must  be  0, 
since  one  point  is  of  Lebesgue  measure  0.  Instead,  consider  a  linear  functional  T  map¬ 
ping  a  smooth  (say  C“)  test  function  m(ar)  having  a  finite  support  (i.e.,  the  domain  where 
it  takes  non-zero  values)  to  the  value  m(0),  i.e.  7^m(z)j=m(0).  This  functional  is  a  well 
defined  entity.  Now,  let  us  agree  to  adopt  a  new  notation  to  express  the  functional; 
write  J 6{x)(  .  )dx,  instead  of  r[  .  ).  As  a  result,  the  above  definition  is  rewritten  as 
j S(x)m(x)dx=m(0).  Thus,  the  delta  function  is  nothing  but  a  notation  for  a  special  func¬ 
tional.  In  fact,  we  do  not  use  the  delta  function  by  itself.  It  is  useful  in  engineering 
problems  only  when  it  is  multiplied  by  some  function  and  integrated.  Hence,  it  suffices 
to  define  only  the  rule  of  integration;  we  need  not  worry  about  its  singularity.  This  is 
the  view  developed  in  detail  by  Schwartz  in  his  theory  of  distributions  [10,  11]. 

TEXTURE  DENSITY  AS  A  FUNCTIONAL 

We  fix  a  window  W  on  the  textured  image  and  define  the  texture  density  J{x,y)  of  a 
dot  texture  over  the  window  W  as  follows: 

Definition  2.1  ^Dot  Density).  The  texture  density  J[x,y)  of  a  dot  texture  over  the  win¬ 
dow  VF  is  a  linear  functional  over  a  set  M,  yet  to  be  specified,  of  test  functions  m(i,j/) 
defined  formally  by 

!j{x,y)m{x.y)dxdy=  J]  H^i.y,),  (o  i) 

PfW  ’ 

where  Pi^Xi,y^  are  thr  dot  texture  elements  on  the  image  plane. 

Since  the  texture  density  is  defined  as  a  functional,  we  need  not  worry  about  the 
singularities  of  f{x,y).  We  can  just  imagine  that  J{x,y)  takes  the  value  infinity  at  texture 
elements  and  0  elsewhere.  All  we  need  is  the  rule  of  integration.  The  texture  density  of 
a  line  segment  texture  is  similarly  defined  as  follows: 

Definition  2.2  (Line  Segment  Density).  The  texture  density  f(x,y)  of  a  line  segment  tex¬ 
ture  in  the  window  IF  is  a  linear  functional  over  a  set  M,  yet  to  be  specified,  of  test 
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functions  m{x,y)  defined  formally  by 

Sj{x,y)m{x,y)dxdy=  Y.  j,  m{x,y)\/ di^+d^,  o) 

l,gw  ‘  '  ■  ' 

where  Li  are  the  line  segments  on  the  image  plane,  and  the  right-hand  side  is  the  sum  of 
line  integrals  along  each  line  segment. 

Here  again,  we  can  imagine  that  J{x,y)  takes  the  value  infinity  along  texture  line 
segments  and  0  elsewhere. 

HOMOGENEITY 

Now,  we  are  in  a  position  to  define  homogeneity  of  the  texture  density.  Let  J{x,y) 
be  the  texture  density  defined  above.  We  would  like  to  say  that  the  texture  is  homo¬ 
geneous  if  J[x,y)!==^c,  but  again,  due  to  singularities,  this  must  be  interpreted  in  the  weak 
sense  or  in  the  sense  of  a  distribution.  Namely,  let  us  agree  that  what  we  mean  by  this 
is  as  follows: 


Definition  2.3  (Homogeneity).  A  texture  density  J{x,y)  is  homogeneous  if 

Jj{x,y)mix,y)dxdy!===icj^m{x,y)dxdy,  (2.3) 

for  test  functions  m{x,y)  of  the  set  M,  yet  to  be  specified,  where  c  is  a  constant  indepen¬ 
dent  of  the  test  functions  m{x,y). 

The  constant  c  can  be  interpreted  as  the  texture  density  in  an  intuitive  sense,  i.e., 
the  “number  of  dots  per  unit  area”  or  the  “length  of  line  segments  per  unit  area”.  If  we 
use  the  definitions  of  eqns  (2.1)  and  (2.2),  our  definition  is  restated  as  follows. 

Lemma  2.1.  If  the  texture  density  is  homogeneous,  we  have  the  following  approxima¬ 
tion  to  integration: 

—  5]  ”*(*112/1)  for  dot  textures 

^  P,iW 

J  m{x,y)dxdy^  2  / _ 

—  5^  j  m{x,y)\  dx^+dy^  for  line  segment  textures 

L,c  w 


Eqn  (2.4)  can  be  viewed  as  the  Monte  Carlo  simulation  of  integration  of  a  test  func¬ 
tion  Tn(x,y],  where  1/c  is  the  “area  per  dot”  or  the  “area  per  unit  length  line  segment”. 
The  interpretation  is  that  the  texture  is  so  homogeneous  that  the  Monte  Carlo  simula¬ 
tion  of  integration  with  respect  to  the  texture  elements  yields  a  good  approximation. 

Remark  2.1.  If  we  choose  as  a  test  function  m{x,y)  the  characteristic  function 


r-y)=^o 


otherwise 


of  a  region  S,  eqn  (2.4)  states  that  the  number  of  dots  or  the  length  of  line  segments  in 
the  region  5  is  approximately  proportional  to  the  area  of  the  region  S  and  that  the 


1-  IS  the  number  or  the  length  of  texture  elements  in  the  region  5  divided  by  its 
ar**^  This  is  the  interpretation  which  most  people  informally  th'.nk  of  as  the  definition 
'^f  homogeneity  Icf.  the  metnod  of  Aloimonos  and  Swain  [8]  recapitulated  in  Section  9). 


Remark  2.2  We  must  mention  here  that  the  definition  of  homogeneity  defined  above 
depends  on  the  choice  of  the  set  M  of  test  functions  m[x,y).  Even  if  the  texture  is  not 
very  dense,  it  can  be  homogeneous  for  very  smooth  test  functions  rn{x,y)  (i.e.,  viewed 
macroscopic  ally).  However,  it  may  not  be  homogeneous  for  rapidly  varying  test  func¬ 
tions  m(x,y)  (i.e.,  viewed  microscopically).  Figuratively  speaking,  we  are  looking  at  the 
singular  texture  density  J{x,y)  through  fitters  m(i,y),  and  the  homogeneity  is  affected  by 
the  ‘coarseness”  of  the  filter  through  which  we  are  looking.  If,  for  example,  we  take 
M={exptir{kx/  a-hly/ b)} ,  assuming  that  the  window  VF  is  a  rectangle  of  size  2aX26,  and 
set  a  certain  threshold  for  the  approximation  of  eqn  (2.3),  we  can  define  the  degree  of 
homogeneity  by  those  [k,I)  satisfying  the  approximation.  However,  we  do  not  go  into  the 
details,  since  what  we  have  described  so  far  is  sufficient  for  the  discussion  to  follow. 


CHANGE  OF  VARIABLES 

Since  the  integration  over  the  texture  is  defined  as  a  functional  by  eqns  (2.1)  and 
(2.2),  we  must  be  careful  when  we  change  the  variables  of  integration.  The  rule  for  the 
usual  integration  does  not  apply  here.  Consider  two  smooth  functions  u{x,y),  v{x,y)  such 
that  the  correspondence  between  (x,y)  and  (u,v)  is  one-to-one,  and  let  2:(u,ti),  y(u,v)  be  the 
inverse.  Suppose  we  use  («,v)  as  new  coordinates.  Let  IT  be  the  domain  on  the  uu-plane 
corresponding  to  the  window  W  on  the  ary-plane.  Define  the  transformed  texture  density 
j{u,v)  also  as  a  functional  by 

f^u,v)m(u,v)dudv=  f^{x,y)m(x,y)dxdy,  (2.6) 

where  function  rn(u,u)  is  defined  by  m(u,v)^m(x(u,v),y{u,v)).  Now,  consider  how  the  new 
density  .^«,v)  acts,  as  a  functional,  on  a  given  test  function  m{u,v). 

First,  consider  a  dot  texture.  Let  points  P,{u,-,u,)  on  the  uu-plane  be  the  images  of 
points  Pi  on  the  ly-plane.  Then, 

fj[x,y)m{x,y)dxdy=  m(i.-,y,)=  m(a:(ui,t;, •),!/(«, ■,!;,))=  ^  »"(«i.v.)-  fo  7^ 

P,iW  P.eW  P,tW  ^  ’ 

This  relation  defines  the  action  of  density  j[u,v),  as  a  functional,  on  the  test  function 
m(u,t;). 

Next,  consider  a  line  segment  texture.  Let  L,  be  the  line  segments  on  the  ui>-plane 
corresponding  to  the  line  segments  L,- on  the  xy-plane.  Then, 

!j{x,y)m{x,y)dxdy=  L  H^<yW dar+difi 

^  L,CW^' 


=  E  X-  m(a< u,  u) , y( u, u))  v/( x^+Vu) du^+2( y„y„) dudv+( z/ -by/) dv^ 

t  /—  ^kr 


L,CW 


^  r  V  V  Vu  «  +2( y„y„) ««-(-( V y„-  v  „ 

=  E  - 7==5 - Vdu‘+dir.  (2.8) 

'  ■  ■  Vu^+v 


L,cw 


dx(u,v)/du,  etc.,  and  ii=du(t)/dt,  v=dv{t)/dt,  where  (u(<),v(<))  is  an  arbitrary 


parameterization  of  individual  line  segments.  Eqn  (2.8)  defines  the  action  of  density 
j{u,v),  as  a  functional,  on  the  test  function  m(u,t;).  Hence,  we  conclude  as  follows. 


Proposition  2.1  (Density  Transformation).  The  transformed  texture  density  /u.v)  is 
formally  given  by 


J[i(u,v),y(u,v))  for  dot  textures 

V  ( Vu^)  u  V2(  x„x„+  y„y„)uv+{  x/+  v' 


J{2{u,v),y(u,v))- 


(2.9) 


for  line  segment  textures 
where  its  action  as  a  functional  is  defined  by  eqns  (2.7)  and  (2.8). 


Remark  2.3.  For  the  usual  integration  of  a  continuous  density,  we  would  have 

Tn{x,y)dxdi/=J^2{u,v),t^u,v))m{x{u,v),y{u,v)) 
so  that  we  would  obtain 


Vu  yv\ 


dudv, 


(2.10) 


K^,v)=Axiu,v),^u,v)) 


Vu 


Xv 

Vv 


(2.11) 


In  sum,  a  continuous  density  is  multiplied  by  the  Jacobian,  which  is  the  magnification 
ratio  of  “area”,  whereas  a  line  segment  density  is  multiplied  by  the  elongation  ratio  of 
“length”,  which  depends  on  the  orientation  of  individual  line  segments,  and  a  dot  den¬ 
sity  is  multiplied  by  the  increase  ratio  of  “number”,  which  is  always  unity,  since  the 
number  of  dots  is  preserved  by  a  continuous  mapping. 


3.  FIRST  FUNDAMENTAL  FORM 

In  this  section,  we  describe  the  3D  shape  of  a  surface  in  the  scene  in  terms  of  the 
image  coordinates  obtained  through  perspective  projection. 

PERSPECTIVE  PROJECTION 


Let  us  fix  a  Cartesian  zj/2-coordinate  system  in  the  scene.  Let  the  if-axis  be  the  opt¬ 
ical  axis  of  the  camera,  and  (0,0,-y),  the  point  on  the  2:-axis  at  distance  /  from  the  xy- 
plar  ?,  be  the  focal  point.  We  adopt  the  camera  model  that  a  point  in  the  scene  is  pro¬ 
jected  to  the  intersection  of  the  ary-plane  with  the  ray  connecting  the  point  and  the  focal 
point  (Fig.  1).  Thus,  the  ly-plane  plays  the  role  of  the  image  plane  and  /  is  the  focal 
length.  It  is  easily  seen  from  Fig.  1  that  the  correspondence  between  the  point  (X,Y,Z) 
in  the  scene  and  the  projected  point  {x,y)  on  the  image  plane  is  given  by 


(3.1) 


SURFACE  DIFFERENTIALS 

Consider  a  smooth  surface  in  the  scene  whose  equation  is  Z=Z(.Y,  F).  This  equation 
coupled  with  eqns  (3.1)  determines  a  one-to-one  correspondence  between  the  points  in 
the  scene  and  the  points  on  the  image  plane  in  the  form  of  ^=^(3:,^),  Y=Y\x,y).  We 
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first  study  how  the  space  coordinates  X,  Y,  Z  change  on  the  surface  by  considering  the 
relationship  among  differentials  dX,  dY,  dZ  taken  along  the  su’-face.  Taking  the 
differentials  of  both  sides  of  eqns  (3.1),  we  obtain 

fdX-xdZ={f+Z)dx,  fd/-ydZ={f+Z)dy.  (3.2) 

Taking  the  differentials  of  both  sides  of  Z=Z{X,Y),  we  obtain 

dZ=PdX+QdY  {P^dZjdX,  Q=dZldY).  (3.3) 

Eqns  (3.2)  and  (3.3)  can  be  viewed  as  a  set  of  simultaneous  linear  equations  in  dX, 
dY,  dZ.  The  solution  is  obtained  in  the  form 

j[J-Px-Qy)  (3-4) 

Here,  all  the  quantities  on  the  right-hand  sides  are  viewed  as  functions  of  the  image 
coordinates  x,  y  through 

Z=Z[X{x,y),Y{x,y)),  P=^{X{x,y),Y{x,y)),  Q=^[X{x,y),Y{x,y)).  (3.5) 

FIRST  FUNDAJvIENTAL  FORM 

Consider  two  points  (x,y),  {x+dx,y+dy)  infinitesimally  far  apart  on  the  image  plane. 
Let  ds  be  the  SD  distance  between  the  corresponding  points  on  the  surface.  Since 
ds~=dX‘^+dY'^+dZ^,  substitution  of  eqns  (3.4)  yields 


Proposition  3.1  (First  Fundamental  Form). 


ds‘^=  S  9ijdx'd3^, 


•j=i 


(3.7) 


where  z*=z,  x^=y  and 


(l-(Px+Qy)//)~ 

9dx,y)=  ,  yzi  i Py)/ /-( FV Q^)xy/f  ^]=g^i(x.y), 

{l-{Px+Qy)//)^ 

(l-{Px+Qy)/J) 


(38) 


Eqn  (3.7)  is  called  the  first  fundamental  form  and  y=(yij)  is  called  the  first  funda¬ 
mental  metric  tensor.  The  first  fundamental  form  of  eqn  ^3.7)  indeed  plays  a  fundamen¬ 
tal  r  j  in  computing  3D  quantities  in  terms  of  the  image  coordinates.  For  example,  con¬ 
sider  an  arbitrary  smooth  curve  L  on  the  image  plane.  The  true  arc  length  of  the 
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corresponding  curve  on  the  surface  is  given  by  integration  /^ds= j=  iffijdx'dx’  on 
the  image  plane. 

Consider  an  infinitesimally  small  square  on  the  image  plane  defined  by  four  points 
(z,t/),  (x+dx,t/),  (x,i/-{-dy),  (x+dx,y+di/).  The  area  of  this  square  on  the  image  plane  is 
dxdy,  but  it  is  easy  to  show  that  the  true  area  of  the  corresponding  region  of  the  surface 
is  given  by  \/det(g}dxdy.  From  eqns  (3.8),  we  obtain 

HPx+Qy)/f  ^  ’ 

Hence,  the  true  area  of  the  region  of  the  surface  corresponding  to  a  region  S  of  the  image 
plane  is  given  by  integration  f^\/det(g)dxdy  on  the  image  plane. 

The  first  fundamental  form  makes  it  possible  to  express  various  other  3D  geometri¬ 
cal  properties  such  as  angle,  Levi-Civita  parallelism,  geodesics  and  surface  gradient  in 
terms  of  the  image  coordinates  (Appendix  A). 

PLANAR  SURFACES 

If  the  surface  is  a  plane  given  by  equation  Z=pX-hqV+r,  where  p,  q,  r  are  con¬ 
stants,  eqns  (3.1)  can  be  solved  for  X,  Tin  the  form 


f-px-qy'  f-px-qy’  f-px-qy 

Since  \+Z/ J={\+r/ f)/[\-{px-\-qy)/ j),  eqns  (3.8)  and  eqn  (3.9)  become  as  follows: 


9uix,y)=-^ — w  X  4  [  ^+p‘--<iy/f+{p‘+ 9^)  r/f  ‘1 , 

{l-{px+qy)/f)* 


(3.10) 


9i2i^,y)=—-] - ^-^■^\-jT[P9M9X+py)/f-{p'^+q-)xy/f'^]=g2^^^^  (3.11) 

{l-(px+qy)/f^^ 


9‘d X, y)=  ;  1 + ?'-2pi//4-(p V q')^/f% 

{\-{px+qy)/f)^ 

\/\+p'+(f[\  +  r/Jl_ 
[\~{px+qy)/ff 


(3.12) 


4.  RECOVERY  OF  PLANAR  SURFACE  ORIENTATION  FROM  TEXTURE 

In  this  section,  we  consider  a  principle  to  compute  the  surface  shape  by  observing 
an  inhomogeneous  texture  density  }{x,y)  on  the  assumption  that  the  true  texture  is 
homogeneous.  We  must  first  study  what  ./(i,y)  looks  like  if  the  true  texture  is  homogene¬ 
ous  Since  the  texture  density  is  defined  as  a  functional,  what  we  need  to  know  is  how 
the  observed  density  ][x,y)  acts  on  a  test  function  m[x,y)  as  a  functional.  Then,  we 
derive  the  basic  equations  to  determine  the  surface  shape  in  term  of  observables  com¬ 
puted  on  the  image  plane. 

DISTORTION  OF  HOMOGENEOUS  TEXTURE 

Consider  temporarily  a  curvilinear  coordinate  system  («, n)  on  the  surface  and 
assume  a  one-to-one  correspondence  n=<i{x.y),  v=v{x.y)  and  r=i(u,n),  y=y{u,v).  Let  IFg 
be  the  region  of  the  surface  corresponding  to  the  window  H’  on  the  image  plane.  Let 
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fQ{u,v)  be  the  homogeneous  texture  density  defined  on  the  surface,  and  let 
mQ{u,v)^Tn{i(u,v),y{u,v)).  According  to  the  assumption  of  homogeneity,  we  have 

Jyyfoi «.  u)  "*o( «,  u)  (iSo^cJ^nioi  u,v)dSo,  (4,1) 

where  dSo  is  the  area  element  of  the  surface. 

Consider  how  eqn  (4.1)  is  expressed  in  terms  of  the  image  coordinates  z,  y.  Since 
the  right-hand  side  is  the  usual  integration  and  dSQ=\/det(g)dxdy,  it  becomes 

^  !^^rMX’yW^^^i9)dxdy.  (4.2) 

The  transformation  of  the  left-hand  side  depends  on  whether  the  texture  consists  of  dots 
or  line  segments. 

For  a  dot  texture,  the  right-hand  side  is  written,  according  to  Proposition  2.1,  as 

P,tW  ’ 

which  can  be  computed  readily  on  the  image  plane. 

For  a  line  segment  texture,  the  right-hand  side  is  written,  in  view  of  Proposition  2.1 
and  ds=y/^^j^ig{jdx'dx‘,  as 

E  L  m{x,y)r{x,y)\/dx-+di/‘, 


where 


r(x.y)=- 


'  9nx  +2gi2iy+922y‘' 


j  .O  .•> 


is  the  elongation  ratio  of  the  line  segment  at  (x,y),  which  depends  on  the  orientation  of 
the  line  segment. 

The  difficulty  is  that  we  cannot  compute  this  integral  on  the  image  plane  unless  we 
know  the  first  fundamental  form,  which  is  dependent  of  the  surface  shape.  Here,  we 
adopt  the  approximation 

n^-y)~(\/3et(^)‘/2. 


The  interpretation  is  that  the  line  segments  are  distributed  nearly  isotropically,  so 
that  if  the  area  is  enlarged  vdet(^  times,  the  individual  line  segments  become  roughly 
(\/det(3))^^"  times  as  long.  Then,  regarding  m(z,y)r(z,j/)  as  a  new  test  function  m{x,y), 
we  can  treat  both  dot  textures  and  line  segment  textures  in  the  same  way.  .Namely,  if 
we  compute,  as  an  observable,  integration 

(-i") 

of  a  test  function  m{x,y)  over  the  window  W,  we  obtain  the  relation 

'^^cj^m{x,y)(\/det(g)Ydxdy,  (4.8) 

where  /c  =  l  for  a  dot  texture  and  /c  =  l/2  for  a  line  segment  texture. 


Remark  4.1.  Eqn  (4  8)  is  interpreted  intuitively  as  follows.  Consider  a  small  region  5 
on  the  image  plane,  and  let  5o  be  the  corresponding  region  on  the  surface.  For  a  dot 


texture,  the  number  of  dots  in  5  is  equal  to  the  number  of  dots  in  5o,  while  the  area  of  5 
is  l/\/det(g)  times  that  of  Sq.  Hence,  the  texture  density  in  S  is  \/det(<?)  times  that  in 
Sq.  For  a  line  segment  texture,  if  the  true  texture  is  nearly  isotropic,  the  total  length  of 
the  line  segments  in  S  is  approximately  l/(\/det((7))'/'  times  that  of  Since  the  area  of 
S  is  l/\/det(^)  times  that  of  5o,  the  texture  density  in  S  is  (\/det(^))^^‘  times  that  in 

GENERAL  PRINCIPLE  OF  SURFACE  RECOVERY 

Our  principle  of  surface  recovery  is  as  follows.  Let  the  object  surface  be  parameter¬ 
ized  so  that  the  procedure  reduces  to  parameter  estimation.  Then,  the  right-hand  side 
of  eqn  (4.8)  is  a  known  form  in  unknown  parameters.  If  we  appropriately  provide  test 
functions  mg(2,!/),  my{x,y),  mo(x,y),  ...,  we  can  compute  the  corresponding  observables  /g, 
/j,  /o,  ...  by  summation  or  integration  on  the  image  plane.  .\s  a  result,  the  necessary 
number  of  equations  are  obtained  in  the  form  of  eqn  (4.8)  to  determine  the  parameter 
values. 

BASIC  EQUATIONS  FOR  PLANAR  SURFACES 

The  simplest  case  is  when  the  surface  is  a  plane.  If  we  replace  the  approximate 
equality  in  eqn  (4.8)  by  equality,  assuming  that  the  true  texture  is  sufficiently  homogene¬ 
ous,  we  obtain  from  eqn  (3.12) 


/  "'^{l-(px+qy)/f}^'^ 


(4.9) 


Now,  provide  three  test  functions  mQ(x,y},  mi(x,y),  vu(x,y),  and  let  7,,  (=0,1,2,  be  the 
corresponding  observables.  If  we  consider  ratios  /j/Zg,  7o/7g,  dropping  off  the  common 
factor  c(\/l-i-p"+q^)''(l^r/7)''',  we  obtain  the  following  equations. 


Proposition  4.1  (Basic  Equations).  The  surface  gradient  (p.q)  is  determined  by  solving 

(=1,2.  (4.10) 


j  m.(  X,  y)-{  /,/  Jq)  mpi  x.  y)  ^ 


{l-{px+qy)/J) 


Zk 


Remark  4.2.  Eqns  (4.10)  are  the  basic  equations  to  determine  (p,q)  and  can  be  solved 
in  principle,  say  by  iterative  search  in  the  p^-space.  Evidently,  three  test  functions  are 
enough  to  determine  the  surface  gradient.  However,  we  can  also  employ  many  more  test 
functions  and  determine  {p,q)  by  some  fitting  scheme  (cf  the  method  of  Dunn  [9]  recast 
in  Section  9). 

SMALL  GRADIENT  APPROXIMATION 

Suppose  the  surface  gradient  (p.q-)  is  close  to  zero  compared  with  the  focal  length  / 
of  the  camera,  so  that  px+qy«f  in  the  window  W.  Since  px+qy=f  is  the  vanishing  line 
or  '‘horizon”  of  the  surface  (vdet(^  of  eqn  (3.12)  becomes  oo  there),  our  assumption  is 
that  the  image  in  the  window  W  is  not  near  the  vanishing  line.  (If  the  vanishing  line 
hnppens  to  be  observed  on  the  image  plane,  its  equation  px-!-7y=/ immediately  tells  us 
the  gradient  (p,q)-)  If  px~^qy«f,  the  Taylor  expansion  around  the  origin  yields 


[\-{px+qy)lff‘ 


-{px+qy)  + 


(4.11) 


If  we  put 


i^^^x{^’y)dxdy,  Ni=  [  yTnlx,y)dxdy 


for  j=0,l,2  and  neglect  higher  order  terms,  the  basic  equations  (4.10)  reduce  to  the  fol¬ 
lowing  linear  equations  in  p,  q: 


^i"('A/'4)M)  '^0)^0  p  / 

M2-{J‘xI  ^2“(-^2/‘A))^0  -1  3/C  ■^o)'^o 


(4.13) 


A  simple  choice  of  the  test  functions  m,{x,y),  /=0,1,2  is 

mo(j:,y)=l,  mi{x:y)=x,  mn{x,y)=y.  (4.14) 

Then,  Jx/ Jq,  Jo/ Jq  nothing  but  the  coordinates  of  the  center  of  gravity  the  texture 
inside  the  window  W.  If  the  window  IT  is  a  rectangle  defined  by  -a<x<a,  -b<y<b, 
the  solution  of  eqns  (4.13)  is  written  as 

p=fx/Ka~,  q=fy/Kb-.  (4.15) 

where  (7,^  is  the  center  of  gravity  of  the  texture  in  the  window  IT.  If  (p,^)=(0,0).  i.e.,  if 
the  surface  is  viewed  orthogonally,  the  center  of  gravity  of  the  texture  coincides  with  the 
origin.  Otherwise,  the  orientation  and  the  magnitude  of  the  “shift”  of  the  center  of 
gravity  of  the  texture  gives  the  surface  gradient  {p,q). 

Remark  4.3  The  accuracy  of  the  result  depends  on  both  the  number  or  length  of  the 
texture  elements  observed  in  the  window  W  and  the  distribution  pattern  of  those  texture 
elements.  Let  .V  be  the  number  or  the  length  of  the  texture  elements  in  the  window  IT. 
The  rule  of  thumb  is  that  the  error  is  approximately  proportional  to  l/\/7v  when  the  tex¬ 
ture  is  completely  random  and  is  approximately  proportional  to  1//Vwhen  the  texture  is 
very  regular  and  periodic  (cf.  Appendix  A).  Textures  we  often  encounter  in  natural 
scenes  and  man-made  objects  are  usually  regular  and  periodic  “tessellations”,  for  which 
high  accuracy  is  expected. 


5.  SCHEME  OF  ITERATIVE  CORRECTION  FOR  PLANAR  SURFACES 

The  method  in  the  previous  section  is  based  on  the  assumption  that  the  gradient 
ip.q)  IS  close  to  zero.  There  exist  methods  which  can  be  applied  when  the  gradient  is  not 
small.  The  first  method  is  based  on  the  following  observation. 

GENERAL  PRINCIPLE 

Suppose  the  camera  is  rotated  by  a  certain  angle  around  its  focus  relative  to  a  sta¬ 
tionary  scene.  .\s  a  result,  a  different  image  is  seen  on  the  image  plane  However,  since 
a  point  on  the  image  plane  actually  corresponds  to  a  “ray”  in  the  3D  scene,  occlusion  is 
not  affected  by  camera  rotation.  If  the  angle  of  camera  rotation  is  known,  the  original 
image  can  be  recovered  as  long  as  the  effect  of  the  image  boundary  is  not  involved.  .\n 
important  fact  is  that  the  image  transformation  due  to  camera  rotation  does  not  require 
any  knowledge  of  the  3D  scene. 


Suppose  the  camera  is  rotated  by  an  orthogonal  matrix  ^=(r„).  The  rotation  of 
the  camera  by  iZ  is  equivalent  to  the  rotation  of  the  scene  by  R~^{=R^.  By  rotation 
R^,  a  point  {X,Y,Z}  moves  to  a  point  {X,Y,2^,  where 


X 

Y 

''ll  ^21  ’’31 

ri2  r22  •'32 

r 

Y 

f+2 

fia  r23  r33 

U+z] 

This  point  is  projected  onto  (i,y)  on  the  image  plane,  where  x=fXj{f+^,  y=fYI[S-\-Z). 
Combining  this  with  eqns  (3.1),  we  obtain  the  transformation  rule 


^  ri3a:+r23t/4-r33/’ 


r  r 22  y+  r^2f 

riiX+r2zy+r2j' 


(5.2) 


Suppose  the  surface  gradient  is  not  small.  We  first  apply  the  method  in  the  previ¬ 
ous  section.  Let  (p,^  be  the  computed  gradient.  This  estimation  may  not  be  accurate. 
Now,  suppose  the  camera  is  rotated  in  such  a  way  that  the  estimated  surface  becomes 
parallel  to  the  image  plane.  As  a  result,  the  gradient  of  the  true  surface  becomes  small, 
so  that  the  method  in  the  previous  section  can  be  applied  again.  Let  {p,q)  be  the  com¬ 
puted  gradient.  If  this  newly  estimated  surface  is  rotated  back  into  the  original  camera 
orientation,  the  obtained  gradient  {p',q')  must  be  a  better  estimate.  This  process  can 
be  applied  repeatedly  until  no  further  improvement  is  obtained.  This  is  the  basic  con¬ 
cept  of  the  iterative  correction  scheme  by  camera  rotation. 

We  should  note  first  that  the  camera  need  not  actually  be  rotated.  Since  the  image 
transformation  due  to  camera  rotation  is  given  by  eqns  (5.2),  the  transformed  image 
f{x,y)  can  be  obtained  by  computation.  However,  we  should  also  note  that  the 
transformed  image  need  not  actually  be  computed.  This  is  because  all  we  need  to  do  is 
apply  appropriate  test  functions  m{x,y)  and  integrate  them  over  the  image.  Since  the 
transformation  is  explicitly  given  by  eqns  (5.2),  the  variables  of  integration  can  be 
changed  so  that  the  integration  can  be  done  over  the  original  image.  In  other  words, 
instead  of  integrating  the  test  functions  m{x,y)  over  the  transformed  image  }[x,y),  we  can 
equivalently  integrate  the  transformed  test  functions  m{x,y)  over  the  original  image 
j{x,y].  Consequently,  our  iterative  correction  scheme  is  performed  by  iteratively  modify¬ 
ing  the  test  functions,  not  the  image. 

CAMERA  ROTATION  SIMULATION 


Let  (p,^  be  the  initial  estimate  of  the  gradient, 
vector  n={ny,n2,n2),  this  means 


In  terms  of  the  surface  unit  normal 


ni=— 


v/l-f-p  ^+q  ^ 
This  vector  makes  angle 


n2==- 


\/l+p  ^+q  ^ 
n=COS"*?l3 


n3= 


1 


\/l+p  ^+'q~ 


(5.3) 

(5.4) 


with  the  unit  vector  lt=(0,0,l)  along  the  z-axis.  The  unit  vector  normal  to  both  n  and  k 
is  given  by 


kXn 


If  the  camera  is  rotated  around  the  vector  I  by  angle  f2  screwwise,  the  estimated 
surface  becomes  parallel  to  the  image  plane.  The  corresponding  rotation  matrix  is  given 
by 

'"n  '‘12 

''12  ^22  ”2  I  (5  6) 

-ni  -no 


where 


—  2  ,  —2 
P  «3+9 

jr2,Tr2  ' 


P‘+?  P 

Hence,  the  image  transformation  is  given  by 

w  .  rnx+ri2y-nj 


'‘12 - _  o  1  o  . 

P  '+9  ‘ 


9  %+P  ~ 
P ’*+9  ‘ 


nii+noy+ns/’ 


yix,y)=f- 


rioX+rn^y-Tt-J 
n^x+rinit'mj  ’ 


and  its  Jacobian  is 


J[x,y)=  -  -  =-/- 

yz  yy 


nii(2;,y)+noi/(x,  !/)-%/■ 

{nxx+n^y+nzff 


The  transformation  of  eqns  (5.8)  maps  the  origin  (0,0)  into  point  If  the  ori¬ 

ginal  image  is  restricted  to  a  window  around  the  origin,  the  transformed  image  is  also 
restricted  to  the  corresponding  domain  IK  around  point  {fp,f^.  Let  {p,q)  be  the  true  gra¬ 
dient  of  the  transformed  surface.  The  Taylor  expansion  around  point  {fp,f^  yields 

7-  '—  ,  .3,c  ^  '  —,3k  i'^^Mi-fP)+^y-f(i))+  ■  ■),  (5.10) 

(l-(px-fgy)/y)3''  (l-pp-qq)^  > 

where 


(5.10) 


j{\-p^qq) 


Q_  3Kq 

A^-phyy) 


(5.11) 


If  eqns  (5.10)  are  substituted,  the  basic  equations  (d.lO)  become 

)q)Mq  Jo)Nq  |-^1 


’  M- 

^i“(4/ 4)'^o 

^2“(4/  4)^0 

where 


4= 


(5.13) 
(5  14) 


Mi=Jj^x-Jp)m,{i.y)dxdy.  Ni=J^^^f^mlx,y)didy.  (5  14) 

and  Jfi.y)  is  the  transformed  texture  density. 

Once  A,  B  are  determined  by  solving  eqn  (5  12),  the  gradient  (p,q)  are  determined 
by  solving  eqns  (5,11),  which  are  rewritten  as 
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■■ 


Ap+Zk/J 

Aq  ■  r 

p 

Bp 

Bq+ZK/f\ 

q  “UJ 

(5.15) 


If  the  computed  gradient  {pt,q)  is  sufficiently  close  to  zero,  the  initial  estimate  is 
correct.  Otherwise,  the  camera  is  rotated  back  into  the  original  orientation,  and  the  sur¬ 
face  gradient  is  transformed  into 


rnP+ruhni 

P  — — : - : - 

riip+Thq+n^ 


ri2P+ro2q-n2 

?  — — : - : - i 

riiP+n^q+Tis 


(5.16) 


These  values  are  better  approximations  to  the  surface  gradient.  This  process  can  be 
iterated  by  setting  p^p  q*—q '  and  repeating  the  previous  procedure  until  convergence. 

COMPUTATION  ON  THE  ORIGINAL  IMAGE 

Consider  how  to  compute  eqns  (5.13)  and  (5.14)  without  actually  transforming  the 
image.  For  a  dot  texture,  we  see,  from  Proposition  2.1,  that  if  we  put 


we  obtain 


mlx,y)=m,{x{x,y),^x,y)),  f=0,l,2. 


P.eW 


(5.17) 

(5.18) 


Thus,  computation  can  be  done  on  the  original  image. 

For  a  line  segment  texture,  we  also  see,  from  Proposition  2.1,  that 


ji=  E  /,  m,{x,y)\/ E{x,y)d^+2F{x,y)dxdy+G[x,y)d]f, 
L,<z  w 


(5.19) 


where 


IT/.  2,-2  /2  l  +  «i^(*(^-y)^+»(a;,y)^-l)-2n,(rii2(x,y)+ri2i;(i,y)) 

— J  - - - 75 - . 

[nxx-^n^y+nzf)^ 

F{x,y)=x^y+yjyy 

_>2  ^l^2(^^.y)^+i<J^.y)^-I)-^i(''i2^(3^.y)+>‘22y(3^,y))-n2(riii(x,y)-hrioy(3:,y)) 

{n,x+n^y+njf  ' 

^—-2,-2  ,2  ^  +  "2^(^^’y)^+i^^>y)^“l)-2"2('’i2^3:,y)-fr22y(a:,y)) 

G[x,y)=Xy  +yy  =/  - - - ^ ^ - - - , 

(Tiii-t-noy-t-nj/)" 

Thus,  computation  can  be  done  on  the  original  image. 

Eqns  (5.14)  are,  on  the  other  hand,  integrations  of  continuous  functions.  Hence,  we 
immediately  obtain 

L{=  x,y)J[x,y)  dxdy, 


Jh/  ^  y)  y)  ^^^y^ 


(5.21) 
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K‘=jJ<^x,y)-f^mlx,y)J{x,y)dxdy. 
These  are  computed  by  a  numerical  integration  scheme. 


Remark  5.1.  Eqn  (5.19)  is  a  rigorous  relation.  A  simple  approximation  is  given, 
corresponding  to  eqn  (4.6),  by 


■4^  E  /,  m,{x,y)\/J{x,y)\/ di'+dr. 


(5.22) 


ALTERNATIVE  METHOD  BY  TAYLOR  EXPANSION 

The  method  described  above  is  somewhat  complicated.  There  exists  an  alternative 
scheme  whose  formulation  is  simpler.  Suppose  (p,5)  is  an  initial  estimate.  We  replace 
eqn  (4.11)  with  the  Taylor  expansion  with  respect  to  {p,q)  at  (p,^: 


(l-(px+9y)/y)3« 
Here,  6p~p-p,  6q=q-'q  and 

L{x,y)= 

Mix,y)= 


^=L{x,y)+M{x,y)6p+N{x,y)6q+ 


{l-{p^+qy)/f}^^ 


Al-{px+qy)/f) 


Mx  v)- 

jiHfx+qyyff'^^  ■ 

Then,  the  basic  equations  (4.10)  become 


■Mi-(yi/yo)A/o  Afi-(/i//o)AQ  rV(A/^)^o' 

_M^(  jy  Jo)Mo  V2-(  J2/ /o)NoJ  Ud  ““  [Lo-{  J2I  Jo)Lo 


where 


(5.23) 


(5.24) 


(5.25) 


^i=jy/nl3:<y)M{x,y)dxdy,  Ni=j^mlx,y)N{x,y)dxdy,  (5.26) 

for  j=0,l,2.  If  dp,  3q  are  sufficiently  close  to  zero,  the  initial  estimate  is  sufficiently  accu¬ 
rate.  Otherwise,  p'=p+6p,  q'=q+8q  are  better  approximations.  This  process  can  be 
iterated  by  setting  p*—p  q*—q  '  and  repeating  this  procedure  until  convergence. 


Remark  5.2.  This  method  is  essentially  the  Newton-Raphson  iteration  of  the  basic 
equations.  Although  this  method  seems  simpler,  the  geometrical  meaning  of  functions 
M{x,y),  N{x,y)  is  not  clear,  while  the  method  of  camera  rotation  has  a  clear  geometrical 
meaning.  In  addition,  the  method  of  Taylor  expansion  still  involves  the  approximation 
(4.6),  while  the  method  of  camera  rotation  is  not  so  much  affected  by  that  approxima¬ 
tion;  the  texture  image  is  exactly  transformed  (cf.  eqns  (5.18),  (5.19))  repeatedly  so  that 
the  true  surface  becomes  more  and  more  parallel  to  the  image  plane,  and  in  the  limit  of 
p— *0,  q—*0,  the  approximation  (4.6)  reduces  to  a  trivial  identity. 


6.  PRINCIPLE  OF  CURVED  SURFACE  RECOVERY  FROM  TEXTURE 

Eqn  (4.13)  can  be  applied  even  if  the  surface  is  not  planar  if  we  choose  windows  of 
an  appropriate  size  so  that  the  surface  is  approximately  planar  in  each  window.  Then, 
eqn  (4.13)  gives  the  gradient  {p,q)  for  each  window.  Thus,  the  global  3D  shape  of  an 
arbitrary  smooth  surface  can  be  recovered  in  principle.  Another  approach  to  curved  sur¬ 
faces  is  to  use  a  surface  model.  Assume,  for  example,  that  the  equation  of  the  surface  is 
given  by 

Z=r+pX+qY+aX  ^+2l3XY+-i  Y'.  (6.1) 

BASIC  EQUATIONS  FOR  CURVED  SURFACES 

If  the  surface  is  planar,  the  observed  texture  in  homogeneity  is  solely  due  to  the  per¬ 
spective  distortion,  whereas  if  the  surface  is  curved,  the  inhomogeneity  is  due  to  two 
separate  sources  -  the  perspective  distortion  and  the  varying  gradient.  In  other  words, 
inhomogeneity  is  not  observed  for  planar  surfaces  if  the  projection  is  orthographic,  i  e  , 
f—*OQ,  while  for  curved  surfaces  inhomogeneity  results  even  if  the  projection  is  ortho¬ 
graphic.  Here,  let  us  consider  orthographic  projection.  Letting  f—*oo,  we  obtain  a^=.V. 
y=  Y  and 

\/  det(  <;)=\/ 1  +p~+q~\/ 1  +AyX+Any+A3X~->-2A^xy+A^i^,  (6.2) 


l+p"+ 


1+P  +?* 


1+p  +<r 


Since  V l+p^+q^  drops  off  when  the  unknown  true  texture  density  c  is  eliminated, 
all  that  can  be  determined  are  parameters  A,-,  »=1,...,5,  which  we  call  texture  density 
parameters.  Consequently,  we  cannot  distinguish  surfaces  which  have  the  same  values 
for  the  texture  density  parameters  A,-,  i=l,...,5. 

If  we  provide  six  test  functions  mo(a:,y),  ...,  m^{x,y)  and  compute  as  observables 

1=0,1, ...,5,  (6.4) 

the  basic  equations  are,  instead  of  eqns  (4.10),  given  by 

Proposition  6.1  (Basic  Equations).  The  texture  density  parameters  are  determined  by 
solving 

ij<Tnlx,y)-{~)mQ{x,y)){y/l+AiX+A2y+A:ir+2A^xy+A^trY<^xdy=0,  i=l,...,5.  (6.5) 


Thus,  six  test  functions  are  enough  to  determine  the  texture  density  parameters  A„ 
j=l,...,5.  Of  course,  we  can  use  many  more  functions  and  determine  the  surface  param¬ 
eters  by  some  fitting  scheme. 
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INTERPRETATION  OF  CURVED  SURFACES 

Once  the  texture  density  parameters  A„  are  obtained,  the  surface  parame¬ 

ters  p,  q,  a,  0,  7  are  determined  by  solving  eqns  (6.3).  The  solution  is  given  as  follows 
(see  Appendix  B  for  proof). 

Proposition  6.2.  The  surface  parameters  p,  q,  a,  0,  7  are  given  by 


P=Ap,  q=kq,  a^ka,  0=k0,  7=A7r. 

(6.6) 

Here,  d,  0,  7  are  given  by 

A3-A5  a 

a=r-l - ,  0= - ,  7=T - , 

16t  8t  16r 

(6.7) 

where 

A2,-\-A^±2y/  A3A^-A4^ 

(6.8) 

Then,  p,  q  are  given  by 

'^Ai~0Ao  .  -0Ai+dtAn 

■4(d7-/?^)  4(01-1-0-) 

(6.9) 

and  k  is  given  by 

ifc-  ^ _ 

^l-p^-q^ 

(6.10) 

AMBIGUITY  OF  INTERPRETATION 

From  eqn  (6.8),  we  see  that  there  exist  at  most  four  solutions.  This  is  an  essential 
characteristic  of  orthographic  projection.  Firstly,  the  projected  image  is  not  affected  if 
we  take  the  mirror  image  of  the  surface  with  respect  to  a  mirror  perpendicular  to  the  z- 
axis  The  four  solutions  consist  of  two  pairs  of  mirror  images.  Another  ambiguity 
occurs  because  the  texture  density  only  tells  about  the  amount  of  surface  inclination  (or 
the  slant)  but  not  the  orientation  of  the  inclination  (or  the  tilt). 

Suppose  the  surface  is  not  a  plane,  i.e.,  parameters  a,  0,  7  are  not  zero  at  the  same 
time.  Still,  there  exist  two  exceptional  cases  (cf.  Appendix  C): 

Case  1.  If  the  surface  has  two  principal  curvatures  of  equal  magnitude,  parameters  a,  0, 
7  are  indeterminate  (i.e.,  parameter  r  of  eqn  (6.8)  becomes  0).  In  this  case,  we  cannot 
tell  whether  the  surface  is  elliptic  (i.e.,  the  Gaussian  curvature  is  positive)  or  hyperbolic 
(i.e.,  the  Gaussian  curvature  is  negative). 

Case  2.  The  four  solutions  for  parameters  d,  0,  7  (two  mirror  image  pairs)  degenerate 
to  two  solutions  (one  mirror  image  pair)  if  and  only  if  the  surface  is  parabolic,  i.e.,  the 
Hessian  a'y-0^  is  zero  (and  hence  the  Gaussian  curvature  is  also  zero).  In  this  case,  the 
inner  square  root  of  eqn  (6.8)  becomes  0. 

Remark  6.1.  Parameters  p,  q  describe  the  gradient  of  the  plane  tangent  to  the  surface 
at  (0,0, r).  They  are  indeterminate  only  in  the  above  two  cases.  In  Case  1,  parameters 
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d,  0,  Tf  are  indeterminate,  so  that  p,  q  are  also  indeterminate.  In  Case  2,  parameters  d, 
0,  7  are  uniquely  determined  except  for  sign.  However,  the  term  a'y-0“  in  the  denomina¬ 
tors  of  eqns  (6.9)  become  0.  Still,  we  can  determine  the  ratio  p\q,  which  indicates  the 
asymptotic  direction  or  the  “ridge”  of  the  surface.  The  fact  that  the  gradient  along  the 
asymptotic  direction  is  indeterminate  is  easily  understood  from  the  orthographic  nature 
of  the  projection;  the  texture  density  is  always  constant  along  this  orientation  all  over 
the  image  plane. 


7.  ALGORITHM  OF  CURVED  SURFACE  RECOVERY  FROM  TEXTURE 

Now,  we  consider  how  to  solve  the  basic  equations  (6.5)  numerically. 

FIRST  ORDER  APPROXIMATION 

If  the  window  W  is  small  and  is  located  at  the  center  of  the  image,  the  Taylor 
expansion  in  x,  y  yields 


{\/T+A^x+A^y+A^rT^lA^iy+A^Y=l->rAx+By-¥C:r+Dxy+E]f+...,  (7.1) 


where 


a=-1a. 


B=—A2, 


C=^[A,+^AY),  D=k{A,+^A,A,),  E=!L[a,+^A,\ 


k-2 


If  we  put 

Mi=J^m,{x,y)dxdy,  Ni=J^mlx,y)dxdy, 

Ri=J^j^m,{x,y)dxdy,  Si=j^ym,{x,y)dxdy,  Ti=J^m,{x,y)dxdy 
for  1=0,1,. ..,5  and  neglect  higher  order  terms,  the  basic  equations  (6.5)  become 
■Mi-(/i/^)Mo  N,-{JJJ^)No  R,-{Ji/Jo)Ro  V(-/i/^)5'o  T,-{J,/Jo)T; 
^2~('4/'A))^H)  ^2~('^2/'^)^0  ^2“('4/'^)^0  ^2~(‘A;/‘^o)'^0  ^2“('4/'^)^0 

^3~( '^3/'A))^0  ■^3"('^3/'A))‘^0  ^3”(■4/'^)^0 

M,-{JJJo)Mo  N,-iJJJo)No  RM/JoWo  S,-(JJJo)S,  T^-[JJ J^)T^. 


LM/Jo)Lo 

Lz-i'^zl  0o)^o 
^4~('A/'4)^0 
'4)^0 


B 

c 

D 

E 


(7.2) 


(7.3) 


(7.4) 


Once  A,  B,  C,  D,  E  are  obtained  by  solving  eqns  (7.4),  the  texture  density  parame¬ 
ters  A„  f=l,2,..,5  are  determined  from  eqns  (7.2)  by 


(7.6) 


Ai - .4, 

AC 


A^=±B, 


2  ^  AC-2  ^2  ^  In  K-2 


K  K‘  /c  AC 


AB,  As=^B-^B^. 


K  AC* 


ITERATIVE  CORRECTION  BY  TAYLOR  EXPANSION 

The  iterative  correction  scheme  by  Taylor  expansion  described  in  the  last  part  of 
Section  5  can  be  applied  to  curved  surfaces  as  well.  Let  A„  be  the  initial  esti¬ 

mates  for  the  texture  density  parameters.  The  Taylor  expansion  at  these  values  yields 

(\/l-l-AiZ-hA2jH-A3ar-t-2A4XjH-A5j^)'‘=mo(ar,y)-l-  m,{x,y)6Ai+  ■  ■  ■  . 

i=l 


(7.6) 


where  M,=A-A„  i=l,...,5,  and 

mo(a:,  !/)=(  i/l  -l-Aii-l-A2iH- A3a:^-l-2A4xy-l- Asj/^)'' , 


mi{x,y)= 


2  [l+AiX+A2y+A3X“ +2A^xy+A^^Y 


m2{x,y)= 


2  (l-t-AiX-f-A2{H-A3r’-f-2A4ZjM-A5tr) 


^3(a;,!/)=-J 


z2 


(7.7) 


2  (l-|-AiZ-l-A2jH-A3Z^-f2A4Z!H-A52/^)‘‘''^^  ’ 


2xy 


2  (l-+-AxZ-l-A2jH-A3ar4-2A4zy-l-A5tr) 


2  [\-\-A^x-\-A2y^A2X^+2A^xy+A^\^y~'^l- 
If  we  use  as  the  test  functions  these  m,(x,y)  themselves  and  put 
•A=/^  X,  y)m,{  X,  y)  dxdy,  i=0, 1 ,  •  •  •  ,5 , 

2/)  y)  » ->=0 , 1 ,  ■  •  ■ ,  5 , 

the  basic  equation  (6.5)  becomes 

I  Aiy  l[M,H-[6i 


(7.8) 

(7.9) 


(7.10) 


where 


(7.11) 


If  all  5A„  1=1,. ..,5,  are  sufficiently  close  to  zero,  the  initial  estimates  are  sufficiently 
accurate.  Otherwise,  A/=A,-t-6A,-  are  better  approximations.  As  before,  this  process  can 
be  iterated  by  setting  A,-^A/  and  repeating  this  procedure  until  convergence. 
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8.  EXAMPLES  OF  COMPUTATION  FOR  SYNTHETIC  IMAGES 

Figs.  2-5  show  texture  images  on  a  planar  surface.  We  take  the  focal  length  /  to 
be  the  unit  of  length.  The  window  size  is  a—b—f  tanl0°=0.176/.  (The  window  is  a  rec¬ 
tangle  of  size  2aX26.)  The  camera  axis  is  assumed  to  pass  through  the  center  of  the 
square  window.  The  true  gradient  is  (p,q')=(  1.500, 0.866)  for  all  the  figures. 

REGULARLY  ALIGNED  DOT  TEXTURE 

Figs.  2a  -  2c  show  projected  images  of  regularly  aligned  dot  textures  on  a  planar 
surface.  If  we  use  mQ{x,y)=l,  mi(x,y)=x,  m^x,y)=y  as  the  test  functions,  i.e.,  if  we 
compute  the  center  of  gravity  of  the  texture,  and  iteratively  correct  these  values  by  the 
method  of  camera  rotation  with  mo(z,y)=l,  m^{x,y)=x-fp,  mo{x,y)=y-fq,  we  successively 
obtain  (p,g)  as  follows. 

Fig.  2a  Fig.  2b  Fig.  2c 

1  (1.459,  1.088)  (1.610,  0.965)  (1.548,  0.958) 

2  (1.402,0.939)  (1.552,0.875)  (1.499,0.871) 

3  (1.397,0.937)  (1.548,0.871)  (1.496,0.867) 

4  (1.397,  0.937)  (1.548,  0.871)  (1.496,  0.867) 

If  we  apply  the  method  of  Taylor  expansion  to  the  same  initial  estimates  with 
mo(i,p)=L(z,i/),  m^{x,y)=M{x,y),  m<i{x,y)=M,x,y),  we  obtain 

Fig.  2a  Fig.  2b  Fig.  2c 

1  (1.459,  1.088)  (1.610,  0.965)  (1.548,  0.958) 

2  (1.549,0.767)  (1.549,0,897)  (1.499,0.869) 

3  (1.527,0.816)  (1.542,0.887)  (1.496,0.864) 

4  (1.529,  0.807)  (1.541,  0.887)  (1.496,  0.864) 

RANDOM  DOT  TEXTURE 

Figs.  3a  -  3c  show  random  dot  textures  on  a  planar  surface.  If  we  compute  the 
center  of  gravity  and  use  the  method  of  camera  rotation,  the  successive  estimates  of  the 
gradient  (p,q)  become 

Fig.  3a  Fig.  3b  Fig.  3c 

1  (1.558,0.798)  (1.662,0.945)  (1.535,0.974) 

2  (1.505,0.695)  (1.617,0.843)  (1.493,0.891) 

3  (1,504,0.694)  (1.613,0.840)  (1.491,  0.888) 

4  (1.504,0.694)  (1.613,0.840)  (1.491,0,888) 

while  the  method  of  Taylor  expansion  yields 

Fig.  3a  Fig.  3b  Fig.  3c 

1  (1.558,0.798)  (1.662,0.945)  (1.535,0,974) 

2  (1.422,0,560)  (1.674,0.832)  (1.527,0.899) 

3  (1.434,0.583)  (1.667,0.834)  (1,525,0.897) 

4  (1.430,0.578)  (1.667,0.834)  (1,525,0.897) 


REGULARLY  ALIGNED  LINE  SEGMENT  TEXTURE 

Figs.  4a  -  4c  show  regularly  aligned  line  segment  textures  on  a  planar  surface.  If  we 
compute  the  center  of  gravity  and  use  the  method  of  camera  rotation,  the  successive 
estimates  of  the  gradient  {p,q)  become 

Fig.  4a  Fig.  4b  Fig.  4c 

1  (1.603,0.994)  (1.762,1.048)  (1.751,1.045) 

2  (1.379,0.747)  (1.534,0.899)  (1.527,0.893) 

3  (1.352,  0.756)  (1.505,  0.875)  (1.500,  0.869) 

4  (1.351,0.758)  (1.504,0.874)  (1.499,0.869) 

while  the  method  of  Taylor  expansion  yields 

Fig.  4a  Fig.  4b  Fig.  4c 

1  (1.603,0.994)  (1.762,1.048)  (1.751,1.045) 

2  (1.569,0.893)  (1.672,0.960)  (1.666,0.958) 

3  (1.568,0.893)  (1.669,0.956)  (1.663,0.954) 

4  (1.568,0.893)  (1.669,0.956)  (1.663,0.954) 

RANDOM  LINE  SEGMENT  TEXTURE 

Figs.  5a  -  5c  show  random  line  segment  textures  on  a  planar  surface.  If  we  com¬ 
pute  the  center  of  gravity  and  use  the  method  of  camera  rotation,  the  successive  esti¬ 
mates  of  the  gradient  {p,q)  become 

Fig.  5a  Fig.  5b  Fig.  5c 

1  (2.821,  0.602)  (2.196,  0.768)  (2.006,  0.910) 

2  (2.275,0.372)  (1.906,0.590)  (1.710,0.722) 

3  (2.088,0.271)  (1.856,0.559)  (1.660,0.684) 

4  (2.111,0.275)  (1.856,0.559)  (1.661,0.684) 

while  the  method  of  Taylor  expansion  yields 

Fig.  5a  Fig.  5b  Fig.  5c 

1  (2  821,0.602)  (2  196,0  768)  (2.006,0.910) 

2  (2.477,0.238)  (2.033,0.561)  (1.825,0.787) 

3  (2.447,0.257)  (2.022,0.567)  (1  820,0.781) 

4  (2.445,0.258)  (2.022,0.567)  (1  820,0.781) 

TEXTURE  ON  A  CURVED  SURFACE 

Fig.  6  is  an  orthographic  view  of  a  regularly  aligned  dot  texture  on  a  quadric  sur¬ 
face.  Here,  we  take  the  window  size  a{=b)  to  be  the  unit  of  length.  The  true  parame¬ 
ters  are  (aa,/?a,7a)=(2,0,0).  (Note  that  parameters  a,  0,  7  have  dimensions  of  1/length.) 
In  this  case,  as  discussed  in  Section  6,  parameters  p,  q  are  indeterminate,  but  parameters 
Q,  3,  7  are  uniquely  determined  except  for  sign.  If  we  use  the  method  of  Section  7  with 
1,  X.  y,  z*,  xy,  y^  as  the  test  functions  for  computation  of  the  initial  estimate  and  apply 
the  method  of  Taylor  expansion,  the  successive  estimates  of  {aa,3b,^a)  become 
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OBSERVATIONS 


i 


Fig.  6 

1  (0.494,  0.000,  -0.007) 

2  (1.066,-0.000,-0.009) 

3  (1.632,  -0.001,  -0.008) 

4  (1  929,-0.001,-0.006) 

5  (1,991,-0.001,-0.006) 


We  see  that  our  methods  can  produce  fairly  good  results.  In  particular,  our  method 
can  be  applied  to  very  sparsely  distributed  textures.  For  very  sparse  textures,  the 
method  of  Taylor  expansion  gives  more  accurate  results  than  the  method  of  camera  rota¬ 
tion.  Otherwise,  both  of  them  yield  almost  the  same  results  for  dot  textures,  but  for  line 
segment  textures  the  method  of  camera  rotation  predicts  more  accurate  results  than  the 
method  of  Taylor  expansion,  as  is  expected  (cf.  Remark  5.2).  The  convergence  is  very 
rapid  for  both  methods.  Only  two  or  three  iterations  are  necessary  to  determine  the  gra¬ 
dient  up  to  two  decimal  places,  and  three  or  four  iterations  up  to  three  decimal  places. 

On  the  whole,  the  results  are  better  for  dot  textures  than  for  line  segment  textures. 
This  is  easily  understood  because  a  line  segment  is  a  coalescence  of  constituent  points  in 
a  restricted  way,  so  that  the  degree  of  homogeneity  is  lower  for  line  segment  textures 
than  for  dot  textures  in  general.  On  the  other  hand,  the  results  are  far  better  for  regular 
textures  than  random  ones,  as  is  also  expected.  This  is  not  a  drawback;  natural  or 
man-made  textures  we  often  encounter  are  usually  “tessellated”  to  a  high  degree  of  regu¬ 
larity.  The  random  texture  shown  here  can  be  regarded,  in  a  sense,  as  the  “worst”  case 
(c*"  Remark  4.3). 


9.  COMPARISON  WITH  OTHER  METHODS 

Let  us  compare  our  method  with  those  of  Aloimonos  and  Swain  [8]  and  Dunn  [9]. 
Both  proposed  schemes  of  surface  shape  recovery  from  texture  based  on  the  geometry  of 
perspective  projection  and  the  assumption  of  texture  homogeneity.  However,  direct 
comparison  is  difficult  because  their  derivations  are  based  on  different  concepts  and 
different  assumptions.  Therefore,  we  now  newly  derive,  in  our  setting,  those  schemes 
which  are  essentially  equivalent  to  (or  actually  better  than)  theirs. 

METHOD  OF  ALOIMONOS  AND  SWAIN 

Consider  three  circular  regions  5o,  5i,  on  the  image  plane  with  centers  (20,2/0)’ 
(zj,2/i),  (xoji/o),  respectively.  Assume  that  they  are  sufficiently  small  compared  with  the 
size  of  the  window  W,  yet  the  texture  is  sufficiently  dense,  so  that  each  region  contains  a 
sufficiently  large  number  of  texture  elements.  Consider,  as  observables,  the  integrals 

/,= J^j{  X,  y)  dxdy,  i=0, 1 ,2 .  ''''.!) 

Since  each  region  5,-  contains  a  large  number  of  texture  elements,  these  integrals  can  be 
approximated,  according  to  eqn  (4.8),  by 

fr«cJ^(\/det((/))''(/a:d2/-  »=0,1,2.  (9.2) 

This  is  equivalent  to  choosing  as  the  test  functions  m,{x,y)  the  characteristic  functions 
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X5,(^.l/)  of  regions  5j  (cf.  Remark  2.1).  On  the  other  hand,  each  region  S,  is  sufficiently 
small,  so  that  the  integral  of  eqn  (9  2)  can  be  replaced  by  the  area  S\  times  the  value  at 
the  center  (x, ■,!/,•),  i.e.. 


/^(\/det((7))'‘didj^5,{\/det(3))''l^^_j,_j^_. 
Thus,  for  a  planar  surface,  the  observables  /,-  are  approximated  by 


.=0,1,2. 

{^-{pxi+qyd/f) 


(9.3) 


(9.4) 


As  a  result,  the  basic  equations  (4.10)  reduce  to 


<(  Hr' m 


fo-Snl  Sit 


^2*^0  1 


-1), 


(9.5) 


from  which  the  gradient  {p,q)  is  determined.  This  is  essentially  the  method  of  Aloimonos 
and  Swain  [8].  They  also  tried  recovery  of  curved  surfaces,  using  a  numerical  relaxation 
scheme. 

Thus,  their  method  requires  that  the  texture  be  locally  homogeneous  in  the  sense 
that  the  texture  is  sufficiently  dense  and  the  homogeneity  condition  is  satisfied  in  arv 
small  region  S  (of  an  appropriate  size).  Furthermore,  a  very  crude  approximation  like 
eqn.  (9.3)  is  used.  In  our  formulation,  however,  the  texture  need  not  be  locally  homo¬ 
geneous.  Some  of  the  textures  shown  in  the  previous  section  are  not  locally  homogene¬ 
ous,  so  that  application  of  the  method  of  Aloimonos  and  Swain  [8]  is  difficult.  For  our 
method,  the  homogeneity  condition  is  required  only  over  the  entire  window  W.  and  the 
integration  is  exactly  performed  over  all  the  texture  elements  in  the  window  W. 

METHOD  OF  DUNN 


Dunn  (9j,  on  the  other  hand,  considered  a  narrow  strip  S  of  width  6  and  length  / 
along  line  x  cos6+y  s\nd=p.  In  our  formulation,  this  process  is  regarded  as  integration  of 
the  texture  density  over  the  strip,  i.e.,  integration  of  the  characteristic  function  \^x,y) 
of  the  strip  5  as  the  test  function  m(x,y).  Take  a  new  jr 'y  ^'coordinate  system  by  rotat¬ 
ing  the  ly-coordinate  system  counterclockwise  by  d  (Fig  7).  Line  x  cosO^y  sm6=p  now 
becomes  x'=p  in  the  new  coordinate  system.  Consider,  as  an  observable,  integration  of 
the  texture  density  over  the  strip  5  (i.e.,  integration  of  the  characteristic  function  \^{i,y) 
of  the  strip  ^  as  the  test  function  m(i,j/)): 

1/2 

K(p.d)=\  J  J{x.y)dx'dy'  (9,6) 

■’  -1/2  M'2 


If  there  exist  a  sufficiently  large  number  of  texture  elements  in  the  strip  .S',  tiie  observ¬ 
able  Kip.B)  IS  approximated,  according  to  eqn  (4.8).  by 

1,2 

K{p.0]^r{  f  ( \/det(  gll'^dx  'dy '  (9  7) 

■'  -i  r  Q-f  -i. 


If  t[|e  approximati'/n  I  f  III  i.s  used  for  a  planar  surface,  this  becomes 


(9.8) 


K(p,0)^cSl(\/l  +p'-^ 1  +  1+'3k{p  cosO^q  smO)p/j). 


.•\s  before,  if  Ki=K{p,,0i),  /=0,1,2,  are  computed  for  three  pairs  of  {p.O)  (i.e.,  for  the 
characteristic  functions  Y5,(^>y)  of  strips  5,  as  the  test  functions  the  basic 

equations  (4.10)  reduce  to 


K  K 

(Picos^i-— PoCOS(9o)p+(pisin^i— — Posin^o)9= 

An 


,  K.  ,  ,  K.  , 

(p^cosOn — — pocos0o)p4-(posin^2 — ^Posin^o)'7=' 


r 

3/C  Aq 

f  K. 

^(1-7^). 
3/C  Aq 


(9.9) 


from  which  p,  q  are  obtained. 

Dunn  |9j,  however,  took  another  approach.  He  searched  for  6  such  that  K{p,9)  does 
not  depend  on  p.  If  6^  is  the  one,  we  see  from  eqn  (9.7)  that  d]  must  satisfy 


p  cos^i+9  sin^i=0, 


(9.10) 


and  hence 

K{pA)^c8I{  i  +  riff^^KY  (9.11) 

IS  a  constant.  Next,  search  for  d  such  that  K(p,6)  has  the  steepest  ascent  with  respect  to 
p  If  Oo  is  the  one,  we  see  from  eqn  (9.7)  that  0o=0i±7r/2  and  that 

p  cosd^+q  sin02='^ P'+T^-  (912) 

and  hence 

K{p.e2)^Ki{l+3KVr+q-p/j).  (9.13) 

The  orientation  of  the  gradient  {p,q)  is  given  by  eqn  (9.9),  and  its  magnitude  is  obtained 
by  computing  the  (average)  gradient  of  /v(p,^2)/^i- 

The  latter  method  can  be  more  robust  to  noise  than  the  former  one,  for  the  former 
method  depends  only  on  three  particular  values  of  K{p,0)  (i.e.,  three  test  functions), 
while  in  the  latter  method  a  large  number  of  strips  (i.e.,  many  different  test  functions,  cf. 
Remark  4.2)  are  used.  Consequently,  some  sort  of  smoothing,  say  fitting  to  a  parametric 
form,  of  h'ip.d)  can  be  used  to  cancel  local  errors  in  the  process  of  searching  for  9i  and 
estimating  the  average  gradient  of  K{p,d.2)l Ky  Dunn  [9l  also  tried  recovery  of  curved 
surfaces  using  a  similar  technique. 

In  any  case,  tne  underlying  approximation  is  essentially  the  small  gradient  approxi¬ 
mation  .-Mthough  the  integration  is  performed  exactly  up  to  linear  approximation  (i.e  , 
eqn  (9  8)),  this  method  also  requires  that  the  texture  be  locally  homogeneous  (though 
mildly  as  compared  with  the  method  of  .Aloimonos  and  Swain  8;),  so  that  the  homo¬ 
geneity  condition  must  be  satisfied  in  each  strip. 


10.  CONCLUDING  REMARKS 

The  methods  proposed  here  have  the  following  salient  features. 

TEXTURE  DENSITY  .\S  .\  FU.NCTION.-VL 
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First,  our  formulation  makes  use  of  the  exact  texture  density  ]{x,y]  having  singulari¬ 
ties,  and  no  smoothing  is  necessary.  This  is  because  we  do  not  use  particular  values  of 
the  texture  density  itself:  all  we  need  is  the  rule  of  integration.  Hence,  the  texture  den¬ 
sity  is  defined  as  a  functional,  or  a  distribution  in  the  sense  of  Schwartz.  This  is  one  of 
the  most  important  differences  from  all  existing  approaches  All  early  methods  assume 
existence  of  a  smooth  texture  density  obtained  by  some  kind  of  local  averaging  and 
directly  use  its  values.  Some  approaches  even  require  the  values  of  the  gradient  of  the 
texture  density  obtained  by  differentiation.  The  use  of  the  values  or  derivatives  of  the 
texture  density  does  not  seem  feasible  in  view  of  the  discrete  nature  of  the  texture. 

DIFFERENTIAL  GEOMETRY  IN  TERMS  OF  IMAGE  COORDINATES 

We  derived  the  exact  relationship  between  the  surface  texture  density  and  the 
observed  texture  density  according  to  the  principles  of  differential  geometry.  The  exist¬ 
ing  methods  seem  to  have  failed  to  obtain  this  exact  relationship.  One  reason,  among 
others,  seems  to  be  that  most  authors  employed  a  certain  intrinsic  “surface  coordinate 
system’’  placed  on  the  surface  as  well  as  an  image  coordinate  system,  trying  to  obtain 
the  rule  of  transformation  from  the  surface  coordinates  to  the  image  coordinates.  This 
usually  results  in  tedious  equations.  The  key  to  success  here  is  the  fact  that  we  do  away 
with  the  surface  coordinates;  all  surface  characteristics  are  described  in  terms  of  the 
image  coordinates  alone,  the  first  fundamental  form  playing  a  fundamental  role. 

COMPUTATIONAL  EFFICIENCY 

Our  method  has  also  an  advantage  from  the  viewpoint  of  computational  efficiency 
The  necessary  data,  or  observables,  are  obtained  by  integration  of  functions  over  the 
image,  and  this  is  essentially  summation  of  the  function  values  over  the  texture  ele¬ 
ments.  Thus,  the  time  complexity  is  simply  <9(/V),  where  iV  is  the  number  of  texture  ele¬ 
ments.  The  access  to  each  texture  element  is  an  independent  process.  This  fact  suggests 
high  speed  performance  by  “parallel  architecture”;  the  image  can  be  divided  in  any  way 
and  the  computation  can  be  performed  independently  and  simultaneously  Although 
Iterations  are  used  in  our  method,  the  convergence  is  very  rapid,  as  was  demonstrated; 
two  or  three  iterations  seem  sufficient. 

PREPROCESSING 

We  should  not  forget  the  fact  that  appropriate  preprocessing  is  necessary,  as  is  also 
the  case  for  any  other  high-level  image  processing.  We  regard  texture  as  composed  of 
dots  without  area  and  line  segments  without  width.  If  the  dots  have  area,  their  centroids 
can  be  used  as  their  positions,  or  their  boundaries  can  be  regarded  as  texture  elements. 
If  the  line  segments  have  width,  their  center  lines  (“skeletons  ’)  or  boundaries  can  be 
regarded  as  texture  elements.  This  is  because  our  method  is  essentially  (weighted) 
number  counting  of  dots  and  (weighted)  length  measuring  of  line  segments.  For 
natural”  texture  images  containing  gray-levels,  a  simple  way  to  do  this  preprocessing  is 
to  just  apply  edge  detection.  Then,  the  detected  “edges”  serve  as  line  segment  texture 
elemen  ts, 

INTEGR.VTIONS  AS  OBSERVABLES 

.Another  advantage  resulting  from  the  use  of  integrations  as  observables  is  that  the 
method  works  for  very  sparse  textures.  .All  existing  methods  including  those  of 
.Aloimonos  and  Swain  8  and  Dunn  9'  have  paid  attention  to  “local”  clues  such  as  the 
number  or  length  of  texture  elements  in  small  regions  of  the  image.  Hence,  the  texture 
must  lie  l-  '-ally  homogeneous:  the  texture  density  must  be  dense  enough  everywhere  so 


that  the  homogeneity  condition  is  satisfied  in  each  of  the  observed  regions.  In  our  for¬ 
mulation,  in  contrast,  the  texture  in  its  entirety  is  observed  directly  through  integration 
over  all  the  texture  elements.  Hence,  the  texture  need  not  be  dense  everywhere,  as  was 
demonstrated  in  the  previous  examples,  and  the  homogeneity  condition  need  be  satisfied 
only  for  the  entire  texture.  (The  idea  of  integrations  as  observables  is  also  used  by 
Kanatani  [12]  for  “shape  from  motion  without  correspondence”.)  On  the  other  hand,  if 
we  gradually  increase  the  texture  density,  our  estimation  approaches  the  numerically 
exact  value.  No  methods  so  far  known  seems  to  have  this  property.  Those  of  Aloimonos 
and  Swain  [8]  and  Dunn  [9]  do  not  have  this  property,  either,  because  various  ad  hoc 
approximations  are  involved. 

DIMENSIONALITY  OF  TEXTURE  ELEMENT 

One  of  the  important  findings  resulting  from  our  analysis  is  the  fact  that  dot  tex¬ 
tures  and  line  segment  textures  cannot  be  treated  in  the  same  manner.  Pixels  constitut¬ 
ing  line  segments  on  the  image  plane  cannot  be  identified  with  pixels  of  a  dot  texture. 
The  necessity  of  this  distinction  does  not  seem  to  have  been  widely  recognized.  All  exist¬ 
ing  methods  including  those  of  Aloimonos  and  Swain  [8]  and  Dunn  [9]  do  not  seem  to 
take  this  effect  correctly  into  account.  It  seems  that  the  texture  density,  whether  of  dots 
or  of  line  segments,  has  been  treated  in  analogy  with  a  continuous  density.  In  this 
paper,  we  have  established  a  rigorous  treatment  of  discrete  densities,  making  a  clear  dis¬ 
tinction  between  dot  textures  and  line  segment  textures.  In  fact,  the  recapitulation  in 
Section  9  is  actually  a  modification  so  that  dimensionality  of  the  texture  elements  is 
correctly  incorporated. 

GENERALITY  OF  THE  PRINCIPLE 

The  main  emphasis  in  this  paper  is  the  generality  of  our  formulation,  from  which 
various  modifications  and  applications  become  possible,  including  the  choice  of  good  test 
functions.  The  methods  of  Aloimonos  and  Swain  [8]  and  Dunn  [9],  for  example,  can  be 
regarded  as  special  variants  of  our  general  principle.  An  important  fact  is  that  our  for¬ 
mulation  can  also  explain  their  methods  and  make  explicit  the  underlying  assumptions 
and  approximations,  while  theirs  can  explain  no  other  methods,  only  their  own.  This 
wide  range  of  flexibility  stems  from  a  mathematically  correct  understanding  of  the 
geometry  of  perspective  projection.  Particular  heuristics  or  ad  hoc  assumptions  and 
approximations  may  result  in  particular  algorithms  which  may  be  useful  sometimes. 
Lacking  generality,  however,  they  usually  do  not  reveal  the  underlying  essential  nature 
of  the  problem  and  are  incapable  of  extension  to  other  problems. 
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APPENDIX  A.  3D  GEOMETRY  IN  TERMS  OF  IMAGE  COORDINATES 

Here  are  some  results  from  differential  geometry  which  are  relevant  in  the  descrip¬ 
tion  of  “shape  from  image”,  describing  3D  properties  in  terms  of  the  image  coordinates. 

INNER  PRODUCT,  NORM,  ANGLE 

Consider  two  vectors  U={U]^,U<^,U^,  V’=(  Vj, Vo,  V3)  tangent  to  the  surface,  making 
angle  d.  Suppose  they  are  “small”  in  the  sense  that  the  first  fundamental  metric  tensor 
jfy  is  almost  constant  along  them.  Let  («*,u^),  (v^,tr)  be  the  projections  of  these  vectors 
onto  the  image  plane;  they  no  longer  make  angle  0.  The  inner  product  of  U  and  V  and 
their  norms  are  computed  from  the  projections  onto  the  image  plane  as  follows; 

( Cl,  V}=  U,  U2 K2+ ^3=  E  (A.  1 ) 

«.;=i 

\\u\\^^/{im,  I|vii=n/(Ty).  (A2) 

Hence,  the  angle  d  is  computed  from  the  projections  of  these  vectors  as 

cose=iU,V)/m\  im|.  (A.3) 

LEVI-CIVITA  PARALLELISM 

Consider  on  the  surface  two  nearby  points  whose  projections  onto  the  image  plane 
are  {x,y],  {x+dx,y+dy].  Again,  consider  two  vectors  of  the  same  length,  “small”  in  the 
sense  described  above,  tangent  to  the  surface  at  these  points.  Let  us  hypothetically  cut 
away  from  the  surface  a  small  patch  which  contains  these  two  tangent  vectors  and  the 
segment  connecting  them  (to  be  precise,  the  developable  defined  as  the  envelope  of  the 
tangent  planes  along  the  segment).  Then,  develop  that  patch,  i.e.,  “  roll  it  out”  on  a 
plane  (Fig.  A).  The  two  tangent  vectors,  having  the  same  length,  are  now  coplanar.  If 
they  are  parallel,  we  say  that  one  of  the  two  tangent  vectors  is  transported  to  the  other 
along  the  segment  parallelly  in  the  sense  of  Levi-Civita.  Let  (u*,u‘),  (u'-t-(5u*,ir-^-(!>ir)  be 
the  projections  of  these  two  tangent  vectors.  The  condition  that  the  corresponding 


(A.4) 


vectors  undergo  the  parallel  transport  in  the  sense  of  Levi-Civita  is  given  by 


where 


is  the  Christoffd  symbol  defined  by 


^9ki  dgjk  v 
dod  ~  daf  ’ 


(A.5) 


and  g~^={g'^)  is  the  inverse  matrix  of  g={gi^.  Assigning  Levi-Civita  parallelism  by  the 
Christoffel  symbol  as  above  is  also  said  to  define  a  Riemann  connection  on  the  image 
plane.  This  enables  one  to  compute  the  curvature  of  the  surface  and  the  curvature  of  a 
curve  on  it  in  terms  of  the  image  coordinates  (the  Gauss -Codazzi  equations  and 
Beltrami’s  formula),  but  we  do  not  go  into  details.  (Refer  to  books  on  differential 
geometry.) 


GEODESICS 


Consider  a  smooth  curve  on  the  surface,  and  let  (2(5), i^s))  be  its  projection  onto  the 
image  plane,  where  the  parameter  s  is  taken  to  be  the  true  arc  length,  defined  by  eqn 
(3.7),  of  the  curve  on  the  surface.  Then,  (i(s),y(s))  is  the  projection  of  the  unit  length 
tangent  vector  to  the  curve.  The  curve  is  called  a  geodesic  if  its  tangent  vectors  are 
always  transported  along  the  curve  parallelly  in  the  sense  of  Levi-Civita.  Hence,  it  fol¬ 
lows  from  eqn  (A.4)  that  the  projection  of  a  geodesic  is  given  by  the  following  differential 
equation: 


i>)+  (A.6) 

It  is  known  that  a  geodesic  is  the  shortest  path  connecting  the  endpoints  along  the  sur¬ 
face. 

GRADIENT 


The  orientation  on  the  surface  Z—Z(X,Y)  along  which  the  value  of  Z  increases  most 
rapidly  on  the  surface  is  indicated  by  vector  (F,Q,P^+Q‘}/(l-hF^+Q‘).  The  projection 
(u\u^)  of  this  vector  onto  the  image  plane  is  given  by 


2 


dP 


(A  7) 


This  vector  indicates  the  orientation  along  which  the  surface  “goes  away”  from  the 
image  plane. 


APPENDIX  B.  ERROR  DUE  TO  RANDOMNESS  OF  THE  TEXTURE 

Consider  a  dot  texture  for  simplicity.  Let  (ii.yi),  ...,  (x/v,y/v)  be  the  coordinates  of 
the  texture  points  in  the  window  -a<i<o,  -6<y<6.  The  center  of  gravity  (x.i/)  is 
givf  .1  by 


First,  consider  the  case  where  the  distribution  is  completely  random.  Suppose  ij, 
are  random  variables  chosen  from  the  uniform  distribution  over  -a<x<a 
independently  from  each  other.  Likewise,  regard  y^,  zs  random  variables  distri¬ 

buted  uniformly  and  independently  over  -b<y<b.  Then,  their  expectation  values  and 
variances  are  given  by 


E[j:,-1=0,  V[J,■)=-^a^  j=1,...,N, 
E[y.l=0,  V|y.)=l62,  ,=l,...,N. 


(B.2) 


It  follows  from  elementary  probability  theory  that  the  expectation  values  and  the  vari¬ 
ances  of  X,  y  are  given  by 


E[zl=0, 

V[J]= 

1  0 

"  3/7“ 

II 

p 

v[y!= 

a 

II 

(B.3) 


Hence,  the  center  of  gravity  is  at  the  origin  on  the  average.  The  magnitude  of  error  is 
estimated  by  the  standard  deviation,  so  that  we  expect  errors  of  about  a/\JZN,  b/\/3N 
for  I,  y,  respectively. 


On  the  other  hand,  consider  another  extreme  case  where  i,-,  i—l,...,N  are  distri¬ 
buted  with  equal  interval  of  distance  2a/N  and  y,-,  »=1 . N  with  equal  interval  of  dis¬ 

tance  2b/ N.  Then,  the  center  of  gravity  must  be  located  within  the  range  of 


a  a 


-±<y<±. 
N  N 


(BA) 


From  eqns  (B.3)  and  (B.4),  we  can  roughly  say  that  the  errors  6x,  by  of  z,  y,  respec¬ 
tively,  are 

(B-5) 

where  l/2<e<l.  The  parameter  e  approaches  1/2  as  the  distribution  becomes  more 
and  more  random,  and  it  c.pproaches  1  as  the  distribution  becomes  more  and  more  regu¬ 
lar. 


APPENDIX  C.  INTERPRETATION  OF  CURVED  SURFACES 

Define  k,  p,  q,  d,  (3,  ^  as  follows: 

h=Vl+p~+q~, 

P=Plk,  q=q/k,  a=a/k,  0=0/ k,  7=7/ A:. 

Eqns  (6.3)  are  now  rewritten  as  follows; 


(C.l) 


Rj«j 


Let  us  define  new  variables 


7 

n 

I 


a=^+i0, 

o’  0 


where  i  is  the  imaginary  unit,  so  that  <t  is  a  complex  number.  Next,  put 

y4-5"i~Ae:  A-l— At  A.A 

T=— - S=— — -+i—,  C.4 

8  8  4  ^  ' 

so  that  S  is  also  a  complex  number.  Then,  the  last  three  equations  of  (C  2)  are 
equivalent  to 

T=i^+(t<x*,  S=2m,  (C.S) 

where  ‘denotes  the  complex  conjugate.  From  the  second  equation,  we  obtain  <T=5/2r. 
Substituting  this  in  the  first  equation,  we  find  that  r  is  the  solution  of 


and  consequently 


or  in  terms  of  A,-,  t=l,...,5, 


t*-Tt^+—SS*=0. 

4 


T^=l(r±\/r‘^-55*), 


r - A3+ A5±2  y/  A3A5-A4*) 


Hence,  eqn  (6.8)  is  obtained,  and  <7  is  given  by  a—S/'lr.  From  eqns  (C.3),  a,  -y  are 
given  by 

a=r+Re[<T],  /^=Im[<r],  7=T--Re[(T],  (C.9) 

from  which  eqns  (6.7)  are  obtained. 

Once  Q,  13,  7  are  obtained,  p,  q  are  determined  from  the  first  two  equations  of  (C  2) 
in  the  form  of  eqns  (6.9).  Finally,  noting 


l+p‘+q'  l+p~+q~  l+p^+q"  kr 

we  obtain  eqn  (6.10).  Thus,  p,  q,  ot,  0,  7  are  given  by  eqns  (6.6). 


APPENDIX  D.  AMBIGUITY  OF  CURVED  SURFACE  RECOVERY 

Consider  the  pathological  cases  ignored  in  Appendix  C.  First,  r  was  assumed  not  to 
be  zero,  since  otherwise  <t=5/2t  is  indeterminate.  From  eqn  (C.7),  t  becomes  zero  if  and 
only  if  5=0  or 


^3—^5' 


Ai=0. 


In  this  case,  eqns  (C.2)  reduce  to 

0(0+^)=O.  (D.2) 

If  r=0,  then  d=(3=']f=0  and  hence  the  surface  is  planar.  Suppose  T^O.  If  j=0,  then 
Q:=±\/T  and  ~i=±\/T'.  If  ;3^0,  then  Q;+Tf=0  and  a~+0‘'=  T  These  two  cases 
correspond  to  surfaces  whose  two  principal  curvatures  have  the  same  magnitude  In 
other  words,  ambiguity  of  d,  0,  7  occurs  for  a  non-planar  surface  if  and  only  if  the  sur¬ 
face  has  principal  curvatures  of  the  same  magnitude,  in  which  case  we  cannot  tell 
whether  the  surface  is  elliptic  (i.e.,  the  Gaussian  curvature  is  positive)  or  hyperbolic  (i.e,, 
the  Gaussian  curvature  is  negative). 

Next,  suppose  we  have  determined  d,  0,  7.  When  solving  the  first  two  equations  of 
(C.2)  for  p,  q,  we  assumed  d7-^^0.  The  condition  that  d7-/^^=0,  or  equivalently 
a'y-^=0,  means  that  the  Hessian  of  the  surface  is  zero  (and  hence  the  Gaussian  curva¬ 
ture  is  also  zero).  This  occurs  for  a  non-planar  surface  if  and  only  if  one  of  the  principal 
curvatures  is  zero  and  hence  the  surface  is  parabolic.  For  the  parabolic  case,  only  the 
ratio  p-.q  is  determined,  indicating  the  asymptotic  direction  of  the  surface.  In  this  case, 
the  proportionality  constant  k  is  indeterminate,  so  that  a,  0,  7  are  also  indeterminate. 
(However,  if  p  and  q  are  known  to  be  small,  we  may  use  the  approximation  of  and 
hence  a^a,  0^0,  7^7.) 

Eqn  (6.8)  indicates  existence  of  four  solutions  They  consist  of  two  pairs  of  mutual 
mirror  images  with  respect  to  a  plane  perpendicular  to  the  z-axis.  The  remaining  ambi¬ 
guity  is  due  to  the  fact  that  the  texture  density  only  tells  about  the  absolute  value  of 
the  slant  and  that  no  information  is  obtained  about  the  ttlt.  From  eqn  (C  7),  this  ambi¬ 
guity  does  not  occur  if  and  only  if  T^=SS*  In  view  of  eqns  (C  5),  this  is  equivalent  to 
T~=(T(T*.  From  eqns  (C.3),  this  is  equivalent  in  turn  to  6''r-i^=0.  In  other  words,  the 
ambiguity  does  not  occur  for  a  non-planar  surface,  except  for  the  mirror  image,  if  and 
only  if  the  surface  is  parabolic. 


(l  *  di,  y  r  dy] 


/  lu-  ^u‘.  (1-  -  fu-'l 


I  -r.  yl  K. 


Fig  A  Two  small  vectors  infinitesimally  far  apart  and  tangent  to  the  surface  are  parallel  in 
the  sense  of  Levi-Civila  if  they  are  parallel  and  of  the  same  length  after  the  small 
patch  containing  them  is  Jeveloped,  i.e.,  “cut  out“  and  “rolled  out”  on  a  plane. 
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